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Computational methods are essential to provide waveforms from coalescing black holes, which 
are expected to produce strong signals for the gravitational wave observatories being developed. 
Although partial simulations of the coalescence have been reported, scientifically useful waveforms 
have so far not been delivered. The goal of the AppleswithApples (AwA) Alliance is to design, 
\^ . coordinate and document standardized code tests for comparing numerical relativity codes. The 

' first round of AwA tests have now being completed and the results are being analyzed. These 

^ initial tests are based upon periodic boundary conditions designed to isolate performance of the 

^SJ . main evolution code. Here we describe and carry out an additional test with periodic boundary 

^ • conditions which deals with an essential feature of the black hole excision problem, namely a non- 

I vanishing shift. The test is a shifted version of the existing AwA gauge wave test. We show how 

, a shift introduces an exponentially growing instability which violates the constraints of a standard 

harmonic formulation of Einstein's equations. We analyze the Cauchy problem in a harmonic gauge 
' and discuss particular options for suppressing instabilities in the gauge wave tests. We implement 

these techniques in a finite difference evolution algorithm and present test results. Although our 
application here is limited to a model problem, the techniques should benefit the simulation of black 
holes using harmonic evolution codes. 
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I. INTRODUCTION 



in 

' Computational methods are essential to provide the waveform from the coalescence of black holes, vi^hich is expected 
^ to produce a strong signal for the gravitational wave observatories being developed. The importance of the binary 
black hole problem to the success of LIGO and LISA has led to major computational efforts, most notably the Binary 
5^1^ _ Black Hole Grand Challenge. Although this Grand Challenge had intermediate successes scientifically useful 

bJ[)' waveforms were not delivered. At present, this remains a problem beyond the reach of any existing code. 

' A recent study |^ of large scale scientific code projects at the Livermore, Los Alamos and Sandia National Labora- 
. tories, funded under the U.S. Department of Energy's Accelerated Strategic Computing Initiative (ASCI), identified 
■ three necessary elements for success: verification, validation and quality management. In the absence of any of those 
^ , three requirements, the report concluded that the results would have little scientific impact because of the impos- 
■ - - ' sibility to judge code reliability. Although the ASCI projects involved highly experienced and qualified teams at 
laboratories with ample resources, only one third of the projects succeeded as planned, another third succeeded later 
than planned and the remaining were eventually abandoned. The failed projects had overly ambitious schedules and 
goals and lacked a conservative methodology that minimized risk. It was expected that the failure rate would have 
been much higher for simulation projects without the experience and resources of the ASCI teams. 

The lessons learned from the ASCI and other studies have been reviewed by Post and Votta j5j. Their experience 
with the peer review process in computational science was that it is not as effective as in experiment or theory, 
because the many possible sources of hidden defects make the referee rely too heavily on plausibility checks as 
opposed to independent reproduction of results. Some of their observations, which are presented below in quotes, 
are very pertinent to numerical relativity: "New methods of verifying and validating complex codes are mandatory 
if computational science is to fulfill its promise for science and society" . The validation of a code implies that the 
predictions of the code are in accord with observed phenomena. For the present status of the binary black hole 
problem, in the absence of any empirical observations, the burden falls completely on verification. 
Post and Votta list five verification techniques, each one having limited effectiveness by itself: 

1. "Comparing code results with an exact answer". 

2. "Establishing that the convergence rate of the truncation error with changing grid spacing is consistent with 
expectations" . 
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3. "Comparing calculated with expected results for a problem especially manufactured to test the code". 

4. "Monitoring conserved quantities and parameters, preservation of symmetry properties and other easily pre- 
dictable outcomes" . 

5. "Benchmarking - that is, comparing results from those with existing codes that can calculate similar problems" . 

The importance of the first four techniques has been recognized by most numerical relativity groups and their 
implementation in practice has improved the integrity of the field. Individual groups cannot easily carry out the 
fifth technique independently. This was the motivation behind formation of the AppleswithApples (AwA) Alliance, 
whose goal is to design, coordinate and document standardized code tests for comparing numerical relativity codes. 
A pivotal step has been the documentation of a first round of AwA tests based upon periodic boundary condition 
(equivalent to a 3-torus with no boundary) [||, designed to isolate performance of the main evolution code. The choice 
of initial AwA tests was biased by considerations of expediency arising from the state of the field at the time. The 
cross fertilization between computational mathematics and numerical relativity was just entering a productive stage. 
Only a few groups had based their codes upon well-posed symmetric or strongly hyperbolic formulations of Einstein's 
equations and fewer groups had an understanding of how to treat boundary conditions. Detailed specifications of a 
second round of tests, involving boundaries (essential for the black hole excision problem), have now been proposed 0. 

Post and Votta emphasized that 

• "Verification and validation establish the credibility of code predictions. Therefore it's very important to have 
a written record of verification and validation results." 

• "Without such programs, computational science will never be credible." 

The first round of AwA tests have now been completed using codes based upon most of the prominent numerical 
relativity formalisms. Results carried out by the participating groups can be viewed in the Alliance CVS depository. 
Instructions for accessing the data in the CVS are available at the Alliance website 0. A second paper is in prepa- 
ration which discusses these first round test results with respect to code performance and improvements in the test 
specifications 

Here we describe and carry out an additional test based upon periodic boundary conditions which deals with an 
essential feature of the black hole excision problem, namely a non- vanishing shift. The test is a shifted version of the 
existing AwA gauge wave test. We detail the test specifications in Sec. ^ 

In Sec. mil we describe some instabilities associated with the gauge wave and shifted gauge wave tests. A discus- 
sion of the Cauchy problem for general relativity depends upon a choice of gauge conditions which reduce Einstein's 
equations to hyperbolic form. See for a recent review. The simplest reduction is in terms of harmonic coordi- 
nates for which well-posedness of the Cauchy problem was first established Previous work 13] has 
revealed that the gauge wave without shift has a constraint preserving instability in harmonic coordinates, i.e the 
gauge wave metric has exponentially growing perturbations which satisfy the harmonic conditions and Einstein's 
equations. We show here how a shift introduces a new type of exponentially growing instability in the standard 
harmonic reduction of Einstein's equations. 

In Sec. lIVI we summarize various options for constructing a hyperbolic evolution algorithm based upon the harmonic 
formulation. Our discussion centers around the standard form of the harmonic formulation adopted in most recent 
analytic treatments. Although the well-posedness of the Cauchy problem guarantees the existence and uniqueness 
of a solution with continuous dependence on the initial data, this is only a necessary condition for computational 
success and is not sufficient. Simulations have shown that instabilities in the gauge wave metrics can quickly introduce 
unacceptable error even with a convergent code |l3l |. In the case of the gauge wave without shift, these instabilities can 
be suppressed by implementing discrete conservation laws obeyed by the principle part of the evolution system |l3l | . 
However, in the case of the shifted gauge wave this technique is not effective by itself because the instability is excited 
by nonlinear terms. In Sec. IIVI we discuss options at the analytic level for suppressing this instability in the shifted 
gauge wave test. One option is to adjust the nonlinear terms in the evolution system by adding terms which vanish 
modulo the harmonic constraints. Another option is the inclusion of harmonic gauge forcing terms, which in principle 
allow the simulation of any nonsingular spacetime region. Gauge forcing terms not only allow the flexibility to "steer" 
around coordinate pathologies that might arise but they also allow the adaptability to carry out standardized tests 
in any specified gauge. 

Section describes the implementation of the harmonic formulation as a finite difference evolution code. We base 
our approach on a code the Abigel code which was developed to implement a well-posed, constraint preserving 
version of the harmonic initial-boundary value problem (IBVP). Here we confine our attention to the Cauchy problem. 
In future work we will extend our test results to the IBVP. Since the preliminary testing of the Abigel code, considerable 
improvement has been made in the underlying numerical techniques. The study of a model nonlinear scalar wave [l5l | 
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shows that semi-discrete conservation laws and summation by parts can be used to formulate stable algorithms. In 
Section [V] we generalize these techniques to the reduced harmonic Einstein equations. These techniques also have 
potential application to the Z4 system jl^ which generalizes the harmonic formulation. 

In Sec. IVII we calibrate the stability, convergence and performance of the code using the gauge wave and shifted 
gauge wave metrics. We show how numerical noise excites instabilities that can be cured by a combination of discrete 
conservation laws and constraint adjustments. 

II. THE SHIFTED GAUGE WAVE TEST 

The standard AwA gauge wave test is based upon the flat metric 

ds^ = (1 - H){-df + dx^) + dy^ + dz^ , (2.1) 

where 

H = H{x-t)=As\ii( ?![(^_^^ (2.2) 



is a sinusoidal wave of amplitude A propagating along the x-axis. In order to test 2-dimensional features, the 
coordinates are rotated according to 

which produces a gauge wave propagating along the diagonal with dependence 

sin — i j' ^ d'^dV2. (2.4) 

Adjusting d or d' to the size of the evolution domain gives periodicity in the x and y directions. 

The shifted gauge wave is obtained from the Minkowski metric ds^ — —dP + dx'^ + djf' + dz^ by the coordinate 
transformation 



(2.5) 
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where d is the size of the evolution domain. This leads to the 4-metric of Kerr-Schild form 

ds^ = ~df + dx^ + dy^ + dz^ + Hkak^dx°'dx>^ (2.6) 

where 



fca = 9a(a;-t) - (-1,1,0,0) (2.7) 

and H is again given by (|2.2|) . This metric describes a shifted gauge wave of amplitude A propagating along the 
X-axis. As above, the coordinate transformation ]2.'6\i rotates the propagation direction to the diagonal. 

The shifted gauge wave test is run in both axis-aligned ID form and diagonal 2D form. We run the test with 
amplitude A = 0.5. We have found that smaller amplitudes are not as efficient for revealing problems. Larger 
amplitudes can trigger gauge pathologies, e.g 5" > (breakdown of the spacelike nature of the Cauchy hypersurfaces), 
more quickly and complicate code comparisons. 

We specify the wavelength c? = 1 in the ID simulations and d' = \f2 in the 2D simulations. We find that at least 50 
grid points lead to reasonable simulations for more than 10 crossing times and therefore make the following choices 
for the computational grid: 

• Simulation domain: 



ID: a; e [-0.5,-^0.5], y = 0, z = 0, d=\ 

2D: a; e [-0.5,-^0.5], y G [-0.5, 4-0.5], z = 0, d' = \/2 
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• Grid: x„ = -0.5 + nda;, n = 0, 1...50p, dx = dy = dz = l/{50p), p= (1,2,4) 

• Time step: dt = dx/i = 0.005//9 

The ID evolution is carried out for T = 1000 crossing times, i.e. 2 x lO^p time steps (or until the code crashes). The 
2D diagonal runs are carried out for T = 100. 

Useful output data are the profiles along the x-axis through the center of the grid (y = z = 0) of gtt, g-xt, and gxx, 
the £2 and £00 norms of the error and the Hamiltonian and momentum constraints, or any other constraints which 
arise in a particular formulation. It is also important to calculate the convergence factor of the numerical error. 



III. INSTABILITY OF THE SHIFTED GAUGE WAVE 



Both the gauge wave metric (|2.1|l and shifted gauge wave metric (|2.6ll are flat vacuum solutions of the Einstein 
equation in harmonic coordinates □x'' = —V^ = 0, where 



= 9"^Kf3 = -^dar^ (3.1) 



and jt'" = y/^gf^". 

Simulation of the AwA gauge wave without shift H2.1II is complicated by the related metric 

dsl = e^*(l - H){-df + dx"^) + dy^ + dz^ (3.2) 

which, for any value of A, is a flat metric which obeys the harmonic coordinate conditions and thus represents a 
harmonic gauge instability of Minkowski space with periodic boundary conditions. 
The shifted gauge wave (|2.6|l has an analogous instability 

dsl = e^\-dt^ + dx^) + dy"^ + dz^ + HKkpdx°'dx^, (3.3) 

which is again a vacuum metric. However, the metric %i.'6\ is not in harmonic coordinates and has a harmonic driving 
term |l3 = f", where 

f " = -XHk". (3.4) 

Thus this particular instability would only be excited in a gauge satisfying H3.4(l . In the standard 3 + 1 description 
x'^ = {t,x^), this determines the propagation of the lapse a = l/\/— 5** and shift /3* — — g**/^" according to 



(3.5) 



and 



t\fh 



where hij is the 3-metric and h — det(/i.y ). 

The exponentially growing metrics H3.2|l and H3.3|l both satisfy the Einstein equations. However, the shifted gauge 
wave also has a different typ e of instability when the evolution system is taken to be the standard harmonic reduction 
of the Einstein tensor 0| 

._ _ y(pp^) l^^M^y^pa^ (3 7) 

which leads to the hyperbolic evolution equation 

E^"" = 0. (3.8) 

(Here is treated formally as a vector in constructing the "covariant" derivative V^F''). In order to see the origin 
of this instability first consider the spatially homogeneous case with amplitude ^ = 0, for which the shifted gauge 
wave metric reduces to the Minkowski metric 77^^/ . Nonlinear perturbations of this metric of the Kerr-Schild form 

a^w = Vi-ii^ + F{t)k^ky (3.9) 
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(where ka = da{x — t) as before) satisfy the reduced harmonic equations H3.8() provided 



E""" = \[{l + F)F,u - FtFt ] k'^k" = 0. (3.10) 



e^*-l. (3.11) 



This has the exponential solution 



The resulting metric 

ffM^ = Vt.u + (e^* - l)k^k^ (3.12) 

solves the reduced equations but unlike H3.3|) it does not solve the Einstein equations because the harmonic condition 
is not satisfied. Instead 

= Ae^'fc'^. (3.13) 

From (|3.7|l . this implies that G*" = 0, i.e. that the Hamiltonian and momentum constraints are satisfied, but that 
the Einstein tensor has the non-vanishing components 

= C ^-^X^e^K (3.14) 

Whether the unstable mode H3.12|l is excited by numerical error depends upon the evolution system. A system 
which enforces G^^ = G^^ = would not excite this mode. In the case of harmonic evolution, the Einstein equations 
are satisfied only indirectly through the satisfaction of the harmonic conditions. These conditions, — 0, are the 
constraints of the harmonic evolution system. Because F^ is not an evolution variable, error of the form H3.13|l can 
be expected to excite the instability. 

The shifted gauge wave (|2.6|) has an exponentially growing instability analogous to (|3.12|) that satisfies the reduced 
harmonic equations but violates the harmonic constraints. The unstable perturbation may be constructed analytically 
by applying the transformation 1)2. 5|l to the exponential solution (|3.12(l of the reduced harmonic equations. Because 
this transformation has the property x — t = x — t, it is simple to verify that F** transforms as a vector. As a result, 
the reduced equations (|3.7|) remain satisfied (since G^''' is a tensor). The resulting metric is 

dsl = ~dt^ + dx'^ + dy^ + (iz^ + ^iJ - 1 + e^*^ kakpdx'^dx'^ , (3.15) 

where now 

Ad f2Tr(x-t)\ 

t.t--cos(^^^j. (3.16) 

The resulting harmonic constraint violation is given by 

F^ = Ae^*V. (3.17) 

The simulation of the shifted gauge wave by a harmonic evolution algorithm based upon the standard reduction H3.7|) 
excites this instability, as exhibited by the results in lVIl However, the instability can be controlled by modifying the 
standard harmonic system, as discussed in the next section. 

IV. THE HARMONIC CAUCHY PROBLEM 

The standard harmonic reduction of the Einstein equations is given by ()3.7|l . with F^ set to zero, which leads to a 
hyperbolic system p. 8(1 which can be cast into the flux conservative form 

+ \g^^'^ (^-±=g'^P{^d^g)dpg + ^gVP^^dpg^P + _L(a^ff)a„5"'3^ = 0. (4.1) 
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The harmonic condition F'^ = comprise the constraints which are sufficient to estabhsh that the Einstein tensor 
vanishes. When the reduced harmonic equations H4.1|l are satisfied, the Bianchi identities imply 

VVaV' + R^T'' ^0, (4.2) 

where the Ricci tensor satisfies 

i?^'' = V^^r"'). (4.3) 

These equations provide the key result that the harmonic conditions propagate in time if the reduced equations are 
satisfied. The historic proof of well-posedness of the initial value problem for Einstein's equations T^l follows from 
the hyperbolicity of 1)3. 7|l and (14.211 . Here hyperbolicity can be defined either in terms of the second differential order 
systems (|4.1f) and (14. 2f) , as in |12|. or by reducing (|4.1|) and (|4.2I) to first order symmetric hyperbolic systems • 
The principal part of E^^'^ , i.e. 

do.{9"^dprn, (4.4) 

has been chosen so that when the remaining nonlinear terms in E'^'^ vanish (or can be neglected) the associated 
conservation laws suppress excitation of the vacuum solution (|3.2I) . These conservation laws, when implemented in 
the discretized system, successfully sup pres s the instability (|3.2() in the simulation of the gauge wave without shift. 
See ^3 for details. The discussion in |T3j | suggests that other fiux conservative forms of the principle part whose 
terms have an analogous conformal weight, such as 

daig^^dfsgf,,), (4.5) 

would be equally effective at suppressing this unstable mode. 

On the other hand, we have found that these conservation laws associated with the principle part H4.4|l are not 
effective in suppressing the instability (|3.15() in the simulation of the shifted gauge wave. In this case, the instability 
is excited by the first derivative terms terms in (|4.1|l which act as a nonlinear source for the principle part (|4.4|l . This 
instability must be handled by a different technique. There are two straightforward generalizations of the standard 
harmonic treatment which leave the principle part of H4.ll) unchanged, so that the well-posedness of the Cauchy 
problem remains intact, but modify the nonlinear terms: (i) the introduction of harmonic gauge source functions or 
(ii) constraint adjustment of the equations. 



A. Harmonic gauge source functions 

Harmonic gauge source functions r''(a;", (7^/3), which are explicit functions of the coordinates and the metric, are 
introduced by replacing the harmonic conditions by = f ^ . The reduced equations then become 

j^t^u ^ ^ y(Mf _ i^/^^Vc^f " = 0. (4.6) 

The generalized harmonic conditions 

:= - f ^ = 0, (4.7) 

are sufficient to ensure that Einstein's equations are satisfied. For this reason we refer to C as the constraints of the 
generalized harmonic formulation. They are related to the standard Hamiltonian and momentum constraints. 
The Bianchi identities imply 

V^VaC" + R^C" = 0. (4.8) 

Because the hyperbolicity of (|4.6|l and (|4.8() is unaffected by such gauge source functions, the Cauchy problem remains 
well-posed. The uniqueness of solutions to (|4.8(l thus ensures that the harmonic constraints C = are satisfied in 
the domain of dependence of the initial Cauchy hypersurface S provided the initial data satisfy C = dtC^ = 0. It is 
straightforward to verify that Cauchy data on S which satisfy the Hamiltonian and momentum constraints = 
and the initial condition C = also satisfy dtC^^ = on 5 by virtue of the generalized reduced harmonic equations 

(ESI). 

The harmonic constraints (|4.7|l also imply equations (|3.5|l and (|3.fci|) governing the time derivatives of the lapse and 
shift. As a result, in addition to the standard Cauchy data, i.e. the intrinsic metric and extrinsic curvature of S 
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subject to the Hamiltonian and momentum constraints, the only other free initial data are the initial choices of lapse 
and shift. 

Thus given Cauchy data 7^" and dt'y^'^ that satisfies the Hamiltonian and momentum constraints in the gauge 
r^' — 0, Cauchy data 7'*" and 9(7'^'' that satisfies the constraints in the gauge = F'' can be obtained by setting 
7^"^ _ ryi^v a,nd dtg^^ = dtg^^ and then determining 9t7*" from 



Finally, 9t7*-' is determined from 



-i=<9^7^'' = -f^. (4.9) 



^ -Af' = dtg'i + ^g'^dtg. (4.10) 



-g' W 



Here the calculation of dig can be expedited using the identity (7**17 = /i, where h — l/det((7*'') is the determinant of 
the 3-metric . Combined with the i-component of (|4.9|) , this yields 



2g h g** 



(f* + -^da'l- (4.11) 



In order to examine whether a gauge source function can be effective in controlling the instability of the shifted 
gauge wave, we consider the simpler problem of spatially homogeneous Kerr-Schild metrics 1)3. 9|) . First consider the 
choice 

f ^ = ct^" (4.12) 
where f^dfj, = dt is the evolution vector. The resulting modification to (|3.1Uf) is 

= i (^(1 + F)F,u - FtFt + cFt^ fc^r = 0. (4.13) 
Because there is no modification to the nonlinear terms, there remain exponentially growing solutions of the form 



F = fce^* - 1 - T- (4-14) 
A 



Another example is the choice 



f^ = ^(7*'^-7fo';), (4.15) 



where we use the notation /[gj = f{0,x^). In the spatially homogeneous case, this choice leads via the constraints 
(lOl to 

5*(7*^-7fo])=-c(7*^-7f,;;), (4.16) 

with solution 7*^ = 7j*^. Thus, if the constraints are satisfied, 7*'^ is forced to retain to its initial value, which has 

the potential advantage of warding off coordinate pathologies. If with the gauge forcing term (|4.15() we make the 
homogeneous Kerr-Schild ansatz (|3.9|) . then the reduced Einstein equations would require 

E''" = ^ ((1 + - F,tF,t + cF[o]F,t^ k^^k^ + ^F* [s'^^l + 5^5"^ = 0. (4.17) 

Again the forcing term has only a linear effect but now it is inconsistent with the Kerr-Schild ansatz (except in the 
trivial case F^t = 0), so that the evolution would in general create other components in the metric and possibly excite 
other instabilities. 

Other gauge source functions can be chosen which introduce nonlinear terms but we have not found an example 
which preserves the Kerr-Schild form 1)3. 9|) . Consequently, an analytic analysis of their effect on the shifted gauge 
wave is difficult to carry out. It is clear from black hole simulations using harmonic coordinates that gauge forcing 
can play a helpful role |23,|21|, but there are few general guidelines to go by. In the case of an analytic testbed, the 
addition of gauge forcing terms beyond those specified in the analytic solution can complicate the test if the resulting 
solution is not known. We have had no computational success in using gauge forcing terms to control instabilities in 
the shifted gauge wave test. 
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B. Constraint adjustment 



There is no extensive knowledge regarding the stabihty of solutions to the constraint systems (|4.2I) or H4.8() . Of special 
importance to numerical evolution is whether constraint violating perturbations of the reduced Einstein equations 
can grow at a fast rate. This question is complicated by the fact that the reduced harmonic equations are not unique. 
They can be adjusted, without affecting their hyperbolicity, according to 



2' 

= £;A"' + v(^f'')-i5^''V„f" + A^'' = (4.18) 

provided the adjustment has the functional form 

^ A^''{x",g"'^,d^g"^,C") (4.19) 
with A^'^ = when C" = 0. The resulting constraint system (|4.8|l becomes 

WWaC + R'iC" - 2V^A^'' = 0, (4.20) 

with the reduced Ricci tensor given by 

jli^u ^ y(p.Qi^) _ + ^S""^- (4.21) 

The standard form of the reduced harmonic equations H3.7|l differ from Fock's Tl'l harmonic formulation, on which 
the original version of the Abigel code was based by the adjustment 

A^" = ^T^'^g'^^^d^g ~ ^g^^T^d^g. (4.22) 
2.9 4g 

In the absence of a general theory, computational experiments are necessary to determine the effect of a given 
adjustment. However, some clues can be provided by the following observations. 
In the linear approximation, with unit lapse and zero shift, the adjustment |22| | 

= -AC'^V^i, A > (4.23) 

leads to the constraint system 

(-92 + v2)C^ = \{^tC^' + S^d^C). (4.24) 



The spatial components satisfy 



(-a^ + V^)(e^*C') = 0, (4.25) 



where r = (e"*"* — 1)/A. For physically reasonable boundary conditions, e.g. periodic boundary conditions, a solution 
F of the wave equation (|4.25|) must have finite energy, so that 

J \VF\^dxdydz < Ki (4.26) 

for some constant Ki. Consequently, 

j \SIC'\^dxdydz < Kie-^^^ (4.27) 

which implies that any initial error in the constraint C" must decay to a spatially constant solution C^(<). Such 
homogeneous solutions of (|4.25|) have the form = K2 + K^e^^^ . By the analogous argument, any constraint 
violation in the time component C* also decays to to a constant value. This result has previously been established for 
the case diC^ ^ using mode analysis p^ . 

Unfortunately, such constraint damping does not extend in a straightforward way to the nonlinear case, where in 
particular it can lead to the excitation of constraint preserving instabilities. For example, this adjustment does not 
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preserve the Kerr-Schild metric type (|3.9|l so that in damping the instabihty H3.15|l it can excite an instabiHty of the 
type 1)3. 2|l . which cannot be controlled by constraint adjustment. See Fig. |^in Sec. I VII for evidence of this behavior. 

The exponential instability (|3.12(l satisfies the reduced harmonic equations but violates the harmonic constraints. 
Although we have found no way to control this instability with harmonic gauge source functions F^, it may be 
suppressed by the constraint adjustment 

A"" = --^C°'dai^g^"'), oO. (4.28) 

The reduced equations (|4.18|) (with F*^ = 0) now give 

2^E^'' = {Fu - FtFt + FFtt + 2cFtFt)k<' k"" = (4.29) 

with solutions 

F=(aii + a2)'/^''^-l. (4.30) 

Thus the exponential growth has been removed. For the case c = 1/2, in which F grows only linearly in time, the 
effect of this adjustment is to remove the lowest differential order terms in from H4.18|l . i.e. 

_ yit^c") + i^'^^V^C" + A^''' = i(-(7''"9„C'' - g^'^dc^C^'). (4.31) 

The growth rate can be completely eliminated by the adjustment 

Af"" ^-'^g'^k'-Tl b>0. (4.32) 

Then H4.18|l gives 

2y/^E'''' = (^{1 + F){F^u + bF^t) - F^tF^t^^k-^ = (4.33) 
with the strongly damped solution 

1 + F = aiexp ^aze"''* j . (4.34) 

The same strongly damped behavior 14.34|l follows from the adjustment 

hC^V f 



where 



°^pa^gpa-\{^pt)'^<yt (4.36) 



g 

is the natural metric of signature (+ + ++) associated with the Cauchy slicing. Here (14.35(1 has the advantage over 
((4.32() that it has a geometric construction which is independent of the Kerr-Schild form of the metric. 

The adjustments (|4.28|) and especially H4.35|l lead to improved performance in the shifted gauge wave tests described 
in Sec IVII Other adjustments can be designed to introduce lowest differential order terms which act as a repulsive 
potential in the wave equation (|4.2UI) . For example, consider the adjustment 

A^"' = -fc.g^'^*^^ fc > 0. (4.37) 

where $ = C"Cq, Then, after contraction with C^, (|4.20|l gives 

Vc«V"$ - 2fc$3 ^ J^, (4.38) 

where T vanishes when W = 0. By choosing k to be large, this repulsive $^ potential might at the least constrain 
C" to be a null vector. 
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V. FINITE DIFFERENCE IMPLEMENTATION 



A first differential order formalism is useful for applying the well developed theory of symmetric hyperbolic systems 
but in a numerical code it introduces extra variables and their associated nonphysical constraints. For this reason, 
we base our code on the natural second order form of the quasilinear wave equations which comprise the reduced 
harmonic system ()4.18|) . They are finite differenced in the flux conservative form 

2^gE^d^{g"Pdpr'')-S^'' = (5.1) 

where S^'^ is comprised of (nonlinear) first derivative terms that do not enter the principle part. 

Numerical evolution is implemented on a spatial grid (cc/, yj, zk) — {Ih, Jh, Kh), < (J, J, K) < N, with uniform 
spacing h, on which a field f{t, x^) is represented by its grid values f[i.j^K]ii) = /(^: UJi zk)- The time integration is 
carried out by the method of lines using a 4th order Runge-Kutta method. We introduce the standard finite difference 
operators Z?oi and D±i according to the examples 

DoxfLJ = -^ifi+ij — fl~l,j) 
D+xfi,j ~ -^{.fi+i,j ^ fi,j) 

the translation operators T±i according to the example 

T±xfi.j = fi±i.j (5.2) 
and the averaging operators A±i and Aoi , according to the examples 

A±xfi.j = l{T±x + l)//,j 



A-oxfi,j — 2^-^+^ T^x)fi,,j- 

Standard centered differences Doi are used to approximate the first derivative terms in (|5.1|l comprising S^'^ . 

We will describe the finite difference techniques for the principle part of 1)5.1(1 in terms of the scalar wave equation 

do.ig^^dp'^) = 0. (5.3) 

Note that the principle part of the linearization of 1)5. 1|) gives rise to (|5.3ll for each component of 7'"'. The non- 
vanishing shift leads to mixed space-time derivatives dtdi which coraplicates the use of standard explicit algorithms 
for the wave equation. This problem has been addressed in [l^ |23L V2M . l25l | . Here we base our work on an evolution 
algorithm shown to be stable for a model 1-dimensional wave equation with shift Provided 5'^ is positive definite, 
as is the case when the shift is subluminal (i.e. when the evolution vector f^da = dt is timelike), this algorithm has 
the summation by parts (SBP) property that gives rise to an energy estimate for (|5.3|) . SBP algorithms have proved 
effective in other numerical relativity codes [IE 113, . 

The algorithm we use is designed to obey semi-discrete versions of two conservation laws obeyed by (|5.3|) . These 
govern the monopole quantity 

g = - / g'"^,^dV. (5.4) 
Jv 

and the energy 

E =]^j^{-g''^l+g'^^,,^,j)dV, (5.5) 

where dV — dxdydz. By assumption, the t — const Cauchy hypersurfaces are spacelike so that —17" > 0. We also 
assume in the following that g^^ is positive definite (subluminal shift) so that E provides a norm. Note that in the 
gravitational case, where <i> represents 7'"^, there are 10 quantities Q°'^ corresponding to 1)5. 4|l . which have monopole, 
dipole or quadrupole transformation properties depending on the choice of indices. 
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For periodic boundary conditions (or in the absence of a boundary), H5.3|) implies strict monopole conservation 
Q t — 0; and, when the cocfhcicnts of the wave operator are frozen in time, i.e. when dtg"^ — 0, H5.3|l impHes strict 
energy conservation E,t — 0. In the time-dependent, boundary-free case, 

Et = lJ^{9f^,c.^,f3)dV. (5.6) 

which readily provides an estimate of the form 

Ej < kE (5.7) 

for some k independent of the initial data for $. Thus the norm is bounded relative to its initial value at t = by 

E < AEoe'''. (5.8) 

The most restrictive value of k depends upon the ratios of the norms of the quadratic forms defined by the integrands 
in H5.5|l and 1)5. 

For the present purpose, it suffices to describe the finite difference evolution algorithm in the 2-dimensional case 
with periodic boundary conditions /o,j — fN,j, fi.o = fi,N- We define the semi-discrete versions of ()5.4|l and 1)5. 5|l as 

N 

Q = h' {-9"'f,t-g''Do,'f). (5.9) 

(/,j)=i 



and 



N 



E = h'^ £, (5.10) 



where 



+ ^{A+ygyy){D+y^)' + ^{A^ygyy){D^y'ff + g-yiDo:.<P)Doy<P. (5.11) 

The energy provides a norm on the discretized system, i.e. E — implies ^ t — D±i^ — (provided the grid spacing 
h is sufficiently small to justify regrouping of terms into perfect squares as in the continuum case). 

The simplest second order approximation to (|5.3(l which reduces in the 1-dimensional case to the SBP algorithm 
presented in is 



where 



W -dt{g''dt<^>) - dt{g'Wo^<^>) - Do,{g'' dt<^>) -Vl<^ = (5.12) 



p2$ ^ iD+,(^iA_,gnD_,^^ +^D^4 iA+,g--)D+,^ 
+ lDJ{A^ygyy)D^y<A+lD^y[ {A+ygyy)D+y'f 



+ Do.l^g'=yDoy^j+Doy\^g'=yDo.<fj. (5.13) 

It follows immediately from the flux conservative form of W that Q,* = for the case of periodic boundary conditions. 
In order to establish the SBP property we consider the frozen coefficient case dtg"^ = 0. Then a straightforward 
calculation gives 

St-^tW = il?+,(^5"$tT„,$t + $tr_,(/'$t) 

+ ^D+J{A_,g--){D_.,^)T_,'i>t) +^D^J (A+,5^-)(i^+,$)r+,$, 
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+ lD + y(^iA^ygyy)iD^y<^>)T^y^t^+^D^y(^iA + ygyy)iD + y^)T+^^^ 

+ ^D+, (^^tT_,{g-yDoy<^) + 5"^(I?o^$)T_,$t 
+ ^D+y (^<^tT-y{g"^Do:,^) + g^y{Do:,^)T^y<^t^ . (5.14) 

Because each term in (|5.14|) is a total D±i, it follows (for periodic boundary conditions) that E^t = when W = 
is satisfied. When the coefBcients of the wave operator are time dependent, an energy estimate can be established 
analogous to H5.8|l for the continuum case. 

We also consider a modification of the algorithm l|5.8|l by introducing extra averaging operators according to 

W ~dt{g"dt<i>) - dt{g''Dor<i>) - D^J {A+,g''){A+,dt<^>)) - = 0, (5.15) 



= {A^^g-^D-.^ + (A+,g^^)D+,<J> 



with 



+ -D + y[{A^ygyy)D^y'^j+-D^y[{A + ygyy)D + y<^ 

+ D^^(^{A+^g''y){A+.,Day<^)^+D^y(^{A+yg''y){A+yDo^-^)y (5.16) 

It is easy to verify that W — W + 0{h?) and both W and W are constructed from the same stencil of gridpoints. 
Although W does not obey the exact SBP property with respect to the energy H5.11|l . the experiments in Sec. IVII 
show that it leads to significantly better performance for the shifted gauge wave test. 

For the time discretization, we apply the method of lines to the large system of ordinary differential equations 

= Ia*.* + -1b*. (5.17) 



obtained from the spatial discretization. Introducing 

h 



= tT, (5.18) 



we obtain the first order system 



t\ _i/baWt 

^)-h\ I j I V 



(5.19) 



We solve the system numerically using a 4th order Runge-Kutta time integrator. 
Dissipation can be added by modifying H5.19|l according to 



Tt ^ Tt + eT/i^I?gI?gT, 

The W algorithm is unstable when g"^^ is not positive-definite or, equivalently, when the evolution direction is 
superluminal (spacelike). There are alternative evolution algorithms for the second differential order ID wave equation 
which are stable in the superluminal case [TEl |2^ . These algorithms have important application to the black hole 
excision problem in treating the region inside the event horizon. However they have no advantage in the shifted 
gauge wave test because the evolution is superluminal only for amplitudes A > \ for which the spacelike nature of 
the Cauchy hypersurfaces breaks down. These alternative algorithms remain stable in the subluminal case but they 
are not as accurate as the W algorithm because they involve a wider stencil of grid points. Nevertheless, it is useful 
to compare their accuracy with the W algorithm. (See Fig. |7|in Sec. IVII ) 

For that purpose we introduce the simplest generalization of the Y algorithm considered in 15] to the 3D case. It 
is related to the W algorithm H5.12|l by 

V = W ^{Vl-Vl))^~Vl,^^ (5.21) 
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where Df^ is defined as in H5.13|l with g'^ replaced by the spatial 3-metric 

h^'^g''-^. (5.22) 
5 

and the shift terms are finite differenced by 



with, for example, 



Equivalently, V and W are related by the 0{h?) terms 



V = W - — [D+.,D.^^-^D+^D^^ + D+yD^y^-^D+yD^y\<^ 



— -Z- Dox.2h n — D+yD-yDoy + D^y 2h 71 — D+xD-xDqx 

2 V 9 9 

+ D+xD^xDox^-^Doy + D+yD^yDoy^-^D^A^. (5.25) 

We also consider the modification 

V = W +{Vl- Vl)^ - Vl2^. (5.26) 

VI. TESTS OF THE EVOLUTION ALGORITHM 

We test the evolution code using the gauge wave and shifted gauge wave testbeds with periodic boundary conditions 
and amplitude A — .5. Test results for the standard AwA gauge wave with amplitudes A = .01 and A — .1 using 
an earlier version of the Abigel code can be found at the Pitt numerical relativity web site . (Those tests were 
performed using an iterated Crank-Nicholson time integrator on a code similar to the W algorithm H5.15|l with the 
constraint adjustments 14.22|l and (|4.28|) with c = .5). The present code shows similar performance for the standard 
AwA gauge wave test. 

The tests are run on grids with N = 50p zones, where p — (1,2,4). We use the norm to measure the error 

f = ||$p-$a„a||oo (6.1) 

in a grid function $p with known analytic value ^ana- We measure the convergence rate at time t 

I |y^4 ~ ^ana \ \ 

using the p = 2 and p = 4 grids (TV — 100 and N — 200). It is also convenient to graph the rescaled error 

£p^^\\^p-^ana\U, (6.3) 

which is normalized to the p = 4 grid. 

A. Tests for the gauge wave without shift 

Figures n and |2] plot our results for the rescaled error £p in g^x, for p = 1,2,4 {N = 50, 100, 200), for the ID and 2D 
gauge wave tests with the gauge wave with amplitude A = .5 using the bare W algorithm (|5.15|) (no additional gauge 
forcing, constraint adjustment or dissipation). The ID runs were stopped at i = 1000. At that time the absolute 
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error for the coarsest grid was unacceptably large (over 100%) but the error for the p = 4 grid was only ~ 20%. We 
found the convergence rates 



"(50) = 2.019, r(500) = 1.677 



(6.4) 



for the error in g^^: rneasured at t=50 and t=500 (corresponding to 50 and 500 grid crossing times). For computational 
economy, the 2D runs, where the wave propagates along a diagonal in the (x, y) plane, were stopped a,t t = 100. We 
measured the convergence rates 



r(10) = 2.084, r(lOO) = 1.568 



(6.5) 



at i = 10 and t = 100. Tests runs with the W algorithm H5.12|l showed similar performance although for the coarsest 
grid the W algorithm had noticeably smaller error, presumably as a result of the extra averaging operators. 
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FIG. 1: Plot of the rescaled error fp(t) for the ID gauge wave. 
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FIG. 2: Plot of the rescaled error Epit) for the 2D gauge wave. 



The bare W or W algorithms give excellent results for these gauge wave tests. No appreciable improvement is 
attained by using the gauge source functions or constraint adjustments discussed in Sec. lIVI or by using dissipation. It 
should be emphasized that this success is due to the discrete conservation laws built into the algorithm. In comparison. 
Fig. O shows snapshots of gtt{x) for a simulation of the gauge wave with A = .5 using an earlier version of the code 
without these conservative properties. The snapshots arc almost perfect matches to the exponentially unstable mode 
llOl . Table n summarizes our results. 
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Test 


Performance 


Non-conservative algorithm, ID 


Exponentially unstable mode (Fig.|^ 


W versus W algorithm, ID 


W has smaller error, due to extra averaging operators 


W versus W algorithm, 2D 


W has smaller error 


W algorithm, convergence, ID 


Long term second order convergence (Fig. 


W algorithm, convergence, 2D 


Long term second order convergence (Fig. ^ 



TABLE I: Summary of tests performed for the gauge wave without shift, with amplitude A — 0.5. Gauge forcing terms, 
constraint adjustments and dissipation all gave no appreciable improvement. 



B. Tests for the Shifted Gauge Wave 

For the shifted gauge wave tests we first consider the bare algorithms (no additional constraint gauge forcing, 
adjustment or dissipation). In this case, the W algorithm shows marked improvement over the W algorithm. This 
is illustrated in Fig. 0] for the ID wave with A=.5, where the W algorithm runs more than twice as long as the W 
algorithm with equal error. For that reason we confine our attention to the W algorithm in the remaining tests. With 
the W algorithm, at t = 50 we measured a convergence rate of 

r(50) = 2.135 (6.6) 

in the error of g^x for the shifted gauge wave. The error is plotted in Fig.Elfor grids with N = 50, 100 and 200. 

The major error in the simulations with the bare algorithms is in the long wavelength exponential mode l|3.12|l . 
which cause the runs with the coarser grids to crash at an early time. This is evident from the series of snapshots 
of gtt{x) shown in Fig. for the N = 100 grid. At i « 200, the long wavelength error triggers the superluminal 
instability of the W algorithm as gu ^ at its peak value. The snapshots of gtt{x) closely match the unstable 
constraint violating mode H3.12|l . The exponential term in (|3.12|) introduces a positive error in gu but a negative 
error in 5" so that the Cauchy hypersurfaces remain spacelike as gu ^0. As a result, this error does not trigger an 
instability in the V algorithm which is stable in the superluminal discussed in Sec. El 

Tests of the bare V and V algorithms (|5.21|) and (|5.26|) with the ID shifted gauge wave gave convergent results 
but the error was considerably larger than with the bare W and W algorithms. Fig. compares the too error in gxx 
between the V and the W algorithms. The V algorithm remains stable well beyond the time t k, 200 when the error 
brings it into the superluminal regime. However, at t ~ 200 the error in the simulation is « 100% and has begun to 
grow rapidly, while the error in the W algorithm is still small. Thus the trade-off for the stability of the V algorithm 
in the superluminal regime is its relative inaccuracy in the subluminal regime compared to the W algorithm. These 
results support the strategy adopted in for a model black hole excision problem, in which the V-algorithm is used 
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FIG. 4: Comparison plots of the unsealed error £{t) in g^^ obtained with the bare W and W algorithms for the ID shifted 
gauge wave test, on an N=200 grid. Although the W algorithm has the SBP property, the additional averaging operators in 
the W algorithm lead to increased accuracy in this nonlinear test. 
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FIG. 5: Plot of the rescaled error £^p(t) in g^^ obtained with the bare W algorithm for the ID shifted gauge wave test. 



in the inner region containing the horizon and the W-algorithm is used in an outer region extending to the outer 
boundary, with the two algorithms blended in a transition region. 

We next consider tests with constraint adjustment. Excellent performance of the W algorithm results from the 
constraint adjustment H4.28|l . The error in g^i^x for the iV = 100 grid is plotted in Fig. Ofor this adjustment with 
c = 0, c = .5 and c = 1. As evident from the plots, the c = 1 run remains under control &t t — 1000 whereas the 
unadjusted run crashes at t w 200. This adjustment with c ~ 0.5 also leads to improvement over the unadjusted case, 
but the results are much more modest. 

The constraint damping adjustment 14.23|l leads to some improvement over the unadjusted case but the results 
are not nearly as effective as the adjustment (|4.28|l . The analysis leading to the decaying behavior H4.25(l in the 
linear regime does not lead to substantially improved performance in this nonlinear test. As shown in Fig. El the 
error obtained with the constraint damped W algorithm, with A = 1, is slightly smaller at early times than for the 
undamped case. However, the error for the damped case then goes through large oscillations and is roughly the same 
as for the undamped case at i = 200. Although the constraint damped case runs longer, the highly oscillatory behavior 
produces unacceptably large error. The constraints exhibit similar oscillation, indicating a coupling to a constraint 
preserving mode. Our experiments indicate significant improvement cannot be obtained by choosing other values of 
the damping coefficient A or by replacing in H4.23() by another timelike vector, such as the vector V"t normal to 
the Cauchy hypersurfaces. (We obtained slightly better results for t".) The addition of numerical dissipation (|5.20|) 
prolongs the run but does not control the large oscillations produced by the constraint damping and does not lead to 
any substantial increase in accuracy. Fig. |51 shows the effect of dissipation with e$ = .1. 

In addition to the c-adjustment (|4.28() . we found that the 5-adjustment H4.35|l also gives excellent performance in 
the shifted gauge wave test. Fig. ^| compares the error, on an iV = 200 grid, obtained using the W algorithm when 
adjusted by (|4.35|l with 6=1 and when adjusted by (|4.28|l with c = 1. Both adjustments are effective in controlling 
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FIG. 6: Snapshots of gtt{x) of the shifted gauge wave obtained with the bare W algorithm on an TV = 100 grid. The behavior 
closely matches the unstable mode H3.12|l . At t ~ 200, the superluminal instability of the W algorithm is triggered. 



W algorithm 
V aleorithm 



FIG. 7: Plots of the unsealed error £{t) in g^x on an A'^ = 200 grid comparing the bare W and V algorithms for the ID shifted 
gauge wave. 



the unstable mode, which is evident in the graph for b = c = 0. The error for the two constraint adjusted cases looks 
very similar, although a closer inspection shows small oscillations in the c— I graph, of the sort evident in Fig. |Slon 
an TV = 100 grid, but which do not appear for b — 1. 

Our experiments with gauge forcing terms, such as (|4.15() . gave no significant improvement in performance. However, 
this is likely a feature of the shifted gauge wave test for which harmonic coordinates do not encounter any pathologies. 
As reported in I2(tII2iI|. we expect gauge forcing terms to be essential for the long term simulation of black holes. 

We found the ID shifted gauge wave tests to be very effective in determining those algorithms which would give 
good performance in the 2D test. For that reason, we limit our presentation of 2D test results for the W algorithm 
with the 6=1 constraint adjustment (|4.35|l . Fig. II II shows the rescaled error obtained using TV — 50, TV = 100 and 
TV — 200 grids. The convergence rate r(lOO) = 2.050 was measured at t = 100. Table UTI summarizes our results. 



VII. CONCLUSION 



Computational fluid dynamics developed into a mature field only after a long struggle during which progress of 
lasting value emerged from the discriminating use of model problems and standardized tests. An important potential 
payoff of numerical relativity is the simulation of gravitational wave sources but sound methodology and testing are 
essential in developing trustworthy codes. The gauge wave tests presented here provide strong caution that a stable 
convergent code is not sufficient to obtain long term simulations, as exhibited by the problems introduced by long 
wavelength unstable modes in Figs. |31and|Sl In the absence of analytic solutions this caution extends to the simulation 
of binary black holes. It would be of value to compare the results reported here for the shifted gauge wave test with 
the performance of codes based upon different hyperbolic reductions of the Einstein equations and different finite 



18 



c=0 

c=0.5 

c=l 



: / 

•■ ; 
: ; 
■■ / 
; 

; / 



: / 
■■/ 



FIG. 8: The unsealed error £{€) in the component g^^, on an = 100 grid, using the constraint adjustment 14.281 in the ID 
shifted gauge wave test. The adjustment with c = \ very eflectively suppresses the instabihty. 




FIG. 9: The effect of constraint damping (A = 1) and dissipation (e = .1) for the ID shifted gauge wave test, on an A'^ = 100 
grid. The graph of the unsealed error E(i) in g:^^ shows a strong oscillation introduced by constraint damping. Although the 
simulation is extended by constraint damping and extended even more by dissipation, the duration of good accuracy is the 
same as with the bare algorithm. The constraints exhibit similar oscillation, indicating coupling of the error to a constraint 
preserving mode. 



c=0, b=0 
c=l, b=0 
c=0, b=l 




FIG. 10: Comparison plots of the unsealed error E{t) in g^^x on an A' = 200 grid for the ID shifted gauge wave test using 
the W algorithm with the constraint adjustment 14.3511 with 6 = 1, the constraint adjustment 14.281 with c = 1 and the bare 
algorithm. The two adjustments show very similar error and both give excellent suppression of the unstable mode excited by 
the bare algorithm. 
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FIG. 11: Plot of the rescaled error £p{t) in g^^ for the 2D shifted gauge wave, with adjustment 6=1. 



Test 


Performance 


W versus W algorithm, ID 


W has smaller error (Fig. ^ 


Stability of bare W algorithm, ID 


Conservation laws do not suppress exponential mode (Fig. |S1 ISl 


W versus V algorithm, ID 


In subluminal case, W has much smaller error (Fig. 


Adjusted W algorithm, with c = 1 (eq.OSl. ID 


Adjustment suppresses the instability (Fig. |5J 


W with A = 1 constraint damping 14.231 and dissipation, ID 


Ineffective, triggers large oscillations (Fig.|2j 


W adjusted with c=l ll428l or 6 = 1 llOsll. ID 


Both adjustments suppress the instability (Fig. IIOII 


W algorithm adjusted with b — 1 (14.3511 . convergence, 2D 


Long term second order convergence (Fig. IllH 



TABLE II: Summary of tests performed for the gauge wave with shift, with amplitude A — 0.5. 




difference approximations. 

Our results show that discrete conservation laws for the principle part of the evolution system are a good starting 
point for designing an algorithm but they do not necessarily control nonlinear instabihties in the analytic problem. 
The same holds true for the standard harmonic reduction of the Einstein tensor. In the case of periodic boundary 
conditions, the techniques introduced in this paper are very effective in suppressing the long wavelength instabilities 
which exist in the shifted gauge wave problem. We have extended these studies of the shifted gauge wave to the 
initial-boundary value problem and a report is in preparation p^ . 
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